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ABSTRACT 

In this paper, we introduce a new approach to constructing unbiased estimators when computing expectations 
of path functionals associated with stochastic differential equations (SDEs). Our randomization idea is 
closely related to multi-level Monte Carlo and provides a simple mechanism for constructing a finite 
variance unbiased estimator with "square root convergence rate" whenever one has available a scheme that 
produces strong error of order greater than 1/2 for the path functional under consideration. 

1 INTRODUCTION 

We have recently developed a general approach to constructing unbiased estimators, given a family of 
biased estimators. It turns out that the conditions guaranteeing its validity are closely related to those 
associated with multi-level Monte Carlo methods; see Rhee and Glynn (2012) for details and a more 
complete discussion of the theory. In this paper, we briefly describe the idea in the setting of computing 
solutions of stochastic differential equations and provide an initial numerical exploration intended to shed 
light on the method's potential effectiveness. As we will see below, the conditions under which our estimator 
produces an algorithm with "square root convergence rate" essentially coincide with the conditions required 
by multi-level Monte Carlo to converge at the same rate. 

In particular, suppose that we wish to compute an expectation of the form a = Ek{X), where X = 
{X{t) : / > 0) is the solution to the SDE 

dX{t) = ix{X{t))dt + G{X{t))dB{t), (1) 

B = iB{t) : ? > 0) is m-dimensional standard Brownian motion, k : C[0,oo) /?, and C[0,oo) is the space 
of continuous functions mapping [0,°°) into R"^. In general, the random variable (rv) k{X) can not be 
simulated exactly, because the underlying infinite-dimensional object X can not be generated exactly. 
Instead, one typically approximates X via a discrete-time approximation Xh{-). For example, the simplest 
such approximation is the Euler time-stepping algorithm given by 

Xh{{k+\)h)=Xh{kh) + ii{Xh{kh))h + a{Xh{kh)){B{k+\)h)-B{kh)) (2) 

that defines X/, at the time points 0,h,2h, with Xf, defined at intermediate values via (for example) linear 
interpolation. Because (2) is only an approximation to the dynamics represented by (1), the rv k{Xh) is only 
an approximation to k{X), and consequently k{Xh) is a biased estimator for the purpose of computing a. 
The traditional means of dealing with this is to intelligently select the step size h and number of independent 
replications R as a function of the computational budget c, so as to maximize the rate of convergence. 
However, as pointed out by Duffie and Glynn (1995), such biased numerical schemes inevitably lead to 
Monte Carlo estimators for a that exhibit slower convergence rates than the "canonical" order c~^/^ rate 
associated with Monte Carlo in the presence of unbiased finite variance estimators. 

However, several years ago, Giles (2008) introduced an intriguing multi-level idea to deal with such 
biased settings that can dramatically improve the rate of convergence and can even, in some settings. 
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achieve the canonical "square root" convergence rate associated with unbiased Monte Carlo. His approach 
does not construct an unbiased estimator, however. Rather, the idea is to construct a family of estimators 
(indexed by the desired error tolerance e) that has controlled bias. In this paper, we show how it is possible, 
in a similar computational setting, to go one step further and to produce (exactly) unbiased estimators. 
The remainder of this paper is organized as follows: We discuss the idea in Section 2 of this paper, while 
Section 3 is devoted to an initial computational exploration of this approach. 

2 THE BASIC IDEA 

We consider here a sequence (X/,,^ : « >= 0) of discrete-time time-stepping approximations to X that are 
all constructed on a common probability space in such a way that: 

i) Ek{Xh„) = Ek{X) + 0{hn) as K ^ 0; 

ii) E\k{Xh„)-k{X)\^ = 0{hl') as hn ^ 

for some r > 0, where 0{f{n)) represents a function which is bounded by some constant multiple of /(•) as 
h„ 0. Assuming, as is often the case for such discretization schemes, that the scheme generates normal 
rv's that are intended to mirror the Brownian increments of the process B driving the SDE (as in the Euler 
scheme (1.2) above), the easiest way to algorithmically obtain an approximating sequence Xh^ to X in 
which the X/j„'s are jointly defined on the same probability space is by successive binary refinement, so that 
hn = 2~". In this setting, the new Brownian motion values (S(y2 ("+^)): j odd) required at discretization 
2-(''+i) can be obtained from the existing values (B{j2-") : 7 > 0) by generating B{{2k+ l)2-("+i)) from 
its conditional distribution given B{k2^") and B{{k+ 1)2^"). On the other hand, one's ability to obtain i 
and ii depends both on the path functional k and on one's choice of discretization scheme. 

In particular, suppose that one has established that the discretization X^ exhibits strong order r. This 
implies that 

Esup{\Xh{kh) -X{kh)\'^'' ■.0<k< [t/h\} = 0{h^'). 

Thus, if k is (for example) a "Lipschitz final value" expectation so that k{x) = g{x{\)) for some Lipschitz 
function g : ii is satisfied. In addition, if k is further assumed to be smooth with \k{X)\ integrable, 

then i is satisfied whenever the discretization X^ is known to be of weak order 1 or higher. It should be 
noted that these conditions are (very) closely related to those that appear in the literature on multi-level 
Monte Carlo for SDEs. 

Note that each of the ^(X2-«)'s is a biased estimator for a = Ek{X). To obtain an unbiased estimator, 
observe that ii) implies the existence of p > such that 

00 

££2"^|^(X2-.)-fc(X2-(„-.))P'-<oo. 
n=l 

Consequently, 

00 

£ 2''P\k{X2-n)-k{X2-(n-i))\^'' < oo a.s. 

from which it follows that 

\k{X2-n)-k(X2-(n-i))\ = 0{2-"P) a.s. 
as n ^ 00, and hence (in view of ii), 

00 

k{X) = k{Xy) + ^ k{X2-n)-k{X2-(n-l)) 

n=l 
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We now introduce a rv A^, independent of B, that takes values in the positive integers and has a 
distribution with unbounded support (so that P{N > n) > for n > 1). For such a rv A'^, 

oo 

Ek{X) = Ek{Xi) + ^ E{k{X2-n)-k{X2-(n-i)))I{N > n)/P{N > n) 

n=l 

N 

= E k{Xi) + £ (^(X2-») - kiX2-in-n))/PiN > n) 



^EZ. 



Note that Z is an unbiased estimator for a. This suggests computing a by generating iid replicates of 
the rv Z. Of course, the "square root" convergence rate of such an estimator is not guaranteed. Given the 
role that finiteness of the variance plays in obtaining such convergence rates, we next study this issue. 

Set A, = k{X2-i) — k{X2-{i-i)) for i > 1 and observe that 

N N N N 

EZ^=Ek'{Xi) +EY,Aj/PiN > if + 2Ek{Xi) Y,MP{N > i)+2E'£ £ ^Aj/{P{N > i)P{N > j)) 

i=\ 1=1 i=\ j=i+l 

= Ek^ (Xi ) + £ EAj/P{N > i) + 2Ek{Xi ) £ A; + 2 £ £ EAiAj/P{N > i) 

!=1 !=1 ;=1 j=i+i 



= Ek^{Xi) + Y,EAj/P{N > i)+2Ek{Xi){k{X)-k{Xi))+2Y,EAi{k{X)-k{X2-i))/P{N > i) 
i=l (=1 

CO 

< Ee (Zi ) + 2Ek{Xi )ik{X)-K^i)) + Y. 0{2-^'') /P{N > i) , 

i=l 

SO that varZ < oo if p(^N > i) ~ c2^^' as / — > oo, for < 7 < 2r (where a, ~ bi means that at/bi — )■ 1 as 
i -)• 00). 

Finally, Glynn and Whitt (1992) prove that "square root convergence rate" ensues if varZ < 00 and 
if the expected computational effort required per replication of Z is finite. The expected computational 
"work" required for each Z is (roughly) given by 

N 

where ti is the incremental effort required to compute k{X2-i) (given k{Xi),. . . ,k{X2-{i-i))), and hence can 
be expressed as 

00 

Y^tiPiN>i). (3) 

1=0 

An approximation to is ti = 2'~^ (the number of additional Gaussian rv's needed to generate X2-i). 
In order that (3) be finite, we require that 7 > 1 . Consequently, a square root convergence rate is ensured 
when 2r > 1 ( in which case we can, for example, choose 7= (1 +2r)/2). 

3 A PRELIMINARY COMPUTATIONAL INVESTIGATION 

In this section, we implement our method and compare it to the multilevel Monte Carlo algorithm suggested 
in Giles (2008). We consider two examples: 
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Geometric Brownian Motion (GBM): The process under consideration here is the solution to 

dX{t) = rX{t)dt + GX{t)dB{t) 
subject to ^(0) = 1, r = 0.05, a = 0.2, the functional k is k{x) =x{\). For this set of parameters, 

=0.104506. 

Cox-Ingersoll-Ross process (CIR): Here, X solves 

dX{t) = K{d -X{t))dt + Gy/xit)dB{t), 

subject to X{0) = 0.04, jc = 5, = 0.04, and cr = 0.25. The functional k is taken here to be k{x) =x(l). 
For this example, 

Ek{X) = 0.04. 

The numerical scheme used to solve each of the above SDE's was the Milstein scheme; see Kloeden 
and Platen (1992). For the above problems, we expect r = 1. For the purpose of this paper, we do not try 
to optimize the distribution of N, and instead choose N so that 

P{N > i) = 2^'/^ 

for j > 1. (In other words, we choose y as the midpoint between 1 and 2r, although any choice in (1,2) 
would provide "square root convergence rate".) 

To compare our method to the multi-level Monte Carlo (MLMC) mehtod, we take the view (as in Giles, 
2008) that the root mean square error (RMSE) e to be achieved by the algorithm is given. Giles (2008) 
provides a complete description of how to construct a MLMC estimator achieving approximate RMSE e; 
we have implemented that version of MLMC here. For our unbiased estimator, we generate independent 
and identically disttibuted (iid) replications of the rv Z until such time as the approximate RMSE is less 
than or equal to e. In other words, our estimator for a is 



1 



N{e) 

Z (4) 



where the Z,'s are iid replicates of Z, and 



N{e)=mf< 




is the first time that the sample RMSE of the sample mean drops below e. We use the stopping rule 
A'^(£) in order to permit easy computation for our estimator, although its use is somewhat unnatural for our 
estimator (since its use induces bias in our estimator). 

For each of our two examples, we provide two tables. The first table for each example concerns our new 
estimator (4); IRE stands for "intended relative error", and k% then corresponds to setting e = (^/100)|a|. 
The 90% confidence interval is then obtained by taking 100 replications of (4) for a given e, computing 
the sample mean and sample standard deviation of the 100 observations, and constructing a confidence 
interval based on the normal approximation. The column corresponding to RMSE is the square root of 
the average, over the 100 observations, of the square of (4) minus EZ. (Thus, RMSE is reporting the 
actual root mean square error of the estimator, rather than the intended RMSE that the estimator has been 
designed to attain asymptotically.) The final column, denoted "work", reports a 90% confidence interval 
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Table 1 : Unbiased estimation for GBM 



IRE 


90% Confidence Interval 


RMSE 


Work 


25% 


0.108318 ± 0.004503 


0.02759 


388.3 ± 23.3 


10% 


0.105360 ± 0.001874 


0.01140 


2850.8 ± 145.7 


5% 


0.105244 ± 0.000916 


0.00560 


11278.3 ± 221.0 


2% 


0.104799 ± 0.000388 


0.00237 


73451.9 ± 1198.8 


1% 


0.104500 ± 0.000174 


0.00105 


291488.1 ± 2157.1 


0.5% 


0.104492 ± 0.000088 


0.00054 


1170663.9 ± 6067.8 



Table 2: MLMC estimation for GBM 



IRE 


90% Confidence Interval 


RMSE 


Work 


25% 


0.103753 ± 0.002693 


0.01635 


3130.3 ± 2.7 


10% 


0.104648 ± 0.001179 


0.00716 


3842.5 ± 10.5 


5% 


0.104908 ± 0.000551 


0.00337 


6378.9 ±31.9 


2% 


0.104465 ± 0.000218 


0.00133 


24901.3 ± 196.0 


1% 


0.104293 ± 0.000129 


0.00081 


93767.1 ± 735.3 


0.5% 


0.104468 ± 0.000060 


0.00037 


377783.8 ± 4248.2 



Table 3: Unbiased estimation for CIR 



IRE 


90% Confidence Interval 


RMSE 


Work 


25% 


0.039914 ± 0.001779 


0.01080 


883.1 ± 61.4 


10% 


0.040272 ± 0.000741 


0.00450 


5548.5 ± 180.7 


5% 


0.040206 ± 0.000342 


0.00208 


22693.3 ± 690.8 


2% 


0.039880 ± 0.000143 


0.00087 


142480.8 ± 1790.4 


1% 


0.039953 ± 0.000069 


0.00042 


563996.7 ± 3609.4 


0.5% 


0.040022 ± 0.000030 


0.00018 


2259485.0 ± 8235.4 



Table 4: MLMC estimation for CIR 



IRE 


90% Confidence Interval 


RMSE 


Work 


25% 


0.040543 ± 0.000973 


0.00593 


3731.2 ± 13.1 


10% 


0.040186 ± 0.000415 


0.00252 


11518.9 ± 44.8 


5% 


0.039880 ± 0.000208 


0.00127 


42503.8 ± 139.4 


2% 


0.039862 ± 0.000100 


0.00062 


265948.6 ± 535.7 


1% 


0.040005 ± 0.000044 


0.00027 


1098672.3 ± 6866.5 


0.5% 


0.040008 ± 0.000023 


0.00014 


4572949.2 ± 5221.9 
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for the expected number of normal rv's generated to construct (4), based on our 100 samples. The second 
table for each example provides a corresponding set of values for the MLMC estimator. 

Our results are reasonably comparable to those associated with MLMC, despite the fact that we have 
done essentially no tuning to optimize the distribution of N. In addition, our estimator is (arguably) easier 
to implement than MLMC, since (in its current form) there are no algorithmic parameters that are estimated 
"on the fly" within the algorithm (in contrast to MLMC). Thus, the unbiased estimators introduced here 
offer a promising computational alternative to MLMC in the presence of SDE numerical schemes having 
a strong order greater than 1/2. 
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